//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>

#include "ephi.hpp"
#include "octree.hpp"

//#define RESOLUTION 0.1
//#define RESOLUTION 0.05
//#define RESOLUTION 0.025
//#define RESOLUTION 0.0125
#define RESOLUTION 0.00625
//#define RESOLUTION 0.003125

class PolywellOctreeBuilder : public OctreeBuilder
{
public:
  PolywellOctreeBuilder (Statics& statics, const vect3d& c1, const vect3d& c2, size_t limitMemoryMB = 100)
    : OctreeBuilder (c1, c2, 0, 2, limitMemoryMB, RESOLUTION), statics(statics)
  {
    addLimit (vect3d(0, 1, 0));
    addLimit (vect3d(1, 0, -1).normal());
    addLimit (vect3d(0, -1, 1).normal());
  }

  class mydifferator : public differator<6>
  {
    public:
      mydifferator(Statics& statics) : differator<6>(vect3d(1e-04, 1e-04, 1e-04)), statics(statics) { }

      void getValues (const vect3d& pos, prec_t *vals)
      {
        vect3d bf, ef;
        statics.getFields (pos, bf, ef);
        vals[0] = bf.x;
        vals[1] = bf.y;
        vals[2] = bf.z;
        vals[3] = ef.x;
        vals[4] = ef.y;
        vals[5] = ef.z;
      }

    private:
      Statics& statics;
  };

  void addValues (const vect3d& pos, CubicScalar *scalars, CubicVector *vectors)
  {
    prec_t f[6];
    prec_t dfdx[6];
    prec_t dfdy[6];
    prec_t dfdz[6];
    prec_t d2fdxdy[6];
    prec_t d2fdxdz[6];
    prec_t d2fdydz[6];
    prec_t d3fdxdydz[6];

#if 0
    memset (dfdx, 0, 6 * sizeof (prec_t));
    memset (dfdy, 0, 6 * sizeof (prec_t));
    memset (dfdz, 0, 6 * sizeof (prec_t));
    memset (d2fdxdy, 0, 6 * sizeof (prec_t));
    memset (d2fdxdz, 0, 6 * sizeof (prec_t));
    memset (d2fdydz, 0, 6 * sizeof (prec_t));
    memset (d3fdxdydz, 0, 6 * sizeof (prec_t));
#endif

    mydifferator diff(statics);
    diff.getValues (pos, f);
    diff.getDifference (pos, dfdx, dfdy, dfdz, d2fdxdy, d2fdxdz, d2fdydz, d3fdxdydz);

    CubicVector bf, ef;

    bf.deserialize (f, dfdx, dfdy, dfdz, d2fdxdy, d2fdxdz, d2fdydz, d3fdxdydz);
    ef.deserialize (f + 3, dfdx + 3, dfdy + 3, dfdz + 3, d2fdxdy + 3, d2fdxdz + 3, d2fdydz + 3, d3fdxdydz + 3);

    vectors[0].add (bf);
    vectors[1].add (ef);
  }

  bool inContact (const vect3d& pos)
  {
    return statics.inContact (pos);
  }

private:
  Statics &statics;
};

static vect3d getv(Octree& tree, const vect3d& pos)
{
  vect3d p = pos;
  const transf3d_pair& t = cube_polywell_transform48 (p);
  vect3d v[2];
  assert(tree.getValues (p, 0, v));
  return t.b * v[0];
}

#define RADIUS 0.15
#define WR 0.035
#define SPACING 0.01

int main (int argc, char *argv[])
{
  Statics statics;
  make_polywell_cube (statics, RADIUS, WR, SPACING, 200000, 2.3e-07);

  PolywellOctreeBuilder b(statics, vect3d(0, 0, 0), vect3d(0.41, 0.41, 0.41));
  Octree tree;
  b.build (tree);
  tree.save ("otree.tree");
  printf ("read = %u\n", tree.load ("otree.tree"));
  printf ("Tree size = %u\n", (unsigned int)tree.size());

  vect3d bf, ef, pos;

  pos = vect3d(0, 0, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

  pos = vect3d(0.099999999, 0, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

  pos = vect3d(0.099999999, 0.0999999999, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

  pos = vect3d(0.05, 0.05, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

  pos = vect3d(0.05, 0.0999999999, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

  pos = vect3d(0.075, 0.075, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

  pos = vect3d(0.025, 0.025, 0);
  pos.print ("pos");
  getv (tree, pos).print ("octree");
  statics.getFields (pos, bf, ef);
  bf.print ("real");

#if 0
  std::string str = tree.getPovrayScript();
  FILE *fp = fopen ("polywell.inc", "w");
  fprintf (fp, "#declare cc = union {\n");
  fwrite (str.c_str(), str.size(), 1, fp);
  str = statics.get_povray_script();
  fwrite (str.c_str(), str.size(), 1, fp);
  fprintf (fp, "}\nobject { cc rotate<0, 360 * clock, 0> }\n");
  fclose (fp);
#endif

  if (1)
  for ( size_t i = 0; i <= 10; i++ )
    {
      {
      Statics rs;
      rs.addStaticElement (new StaticOctree48 (tree, true, true, false));
      ElectroDynamics dyn(rs);
      Scene scene(dyn);
      Screen screen(800, 800, 0.8, 0.8);
      screen.set_axis (Screen::PLANE_XY, i * 0.05 - 5 * 0.05);
      scene.calc (screen);
      scene.set_coloring (32, true);
      scene.render_map (screen);
      char fn[128];
      sprintf (fn, "otree_bf_%02u.bmp", (unsigned int)i);
      screen.write (fn);
      scene.set_coloring (0, true);
      scene.render_lic (screen);
      sprintf (fn, "otree_bf_lic_%02u.bmp", (unsigned int)i);
      screen.write (fn);
      }
      {
      Statics rs;
      rs.addStaticElement (new StaticOctree48 (tree, true, true, false));
      ElectroDynamics dyn(rs);
      Scene scene(dyn);
      Screen screen(800, 800, 0.8, 0.8);
      screen.set_axis (Screen::PLANE_XY, i * 0.05 - 5 * 0.05);
      scene.calc (screen, Scene::CALC_EFIELD);
      scene.set_coloring (32, true);
      scene.render_map (screen);
      char fn[128];
      sprintf (fn, "otree_ef_%02u.bmp", (unsigned int)i);
      screen.write (fn);
      scene.set_coloring (0, true);
      scene.render_lic (screen);
      sprintf (fn, "otree_ef_lic_%02u.bmp", (unsigned int)i);
      screen.write (fn);
      }
    }

  {
    Statics statics;
    statics.addStaticElement (new StaticOctree48 (tree, true, true, false));
    ElectroDynamics dyn(statics);
    dyn.inject_electron (vect3d(-0.3, 0.003, 0), vect3d(1, 0, 0), 0);
    dyn.setDTM (0.1);
    Screen screen(800, 800, 0.8, 0.8);
    screen.draw_axes();
    statics.draw (screen);
    for ( size_t i = 0; i < 200000; i++ )
      {
        dyn.step();
        screen.set_pixel (dyn.get_pos(0), Screen::RED);
      }
    screen.write ("otree_epath.bmp");
  }

  return 0;
}

